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STORM  HYDROGRAPHS  BY  TWO -STAGE  CONVOLUTION 
DETAILED  FORM  OF  COMPUTATIONS 


W.  M.  Snyder 


Convolution  is  a  conventional  method  of  input-output  analysis  in 
hydrology  (1) .  Two-stage  convolution  has  been  used  in  hydrologic  models 
in  the  Southeast  Watershed  Laboratory  (2,  3,  4,  5).  The  method  is  useful 
for  both  prediction  and  identification  as  discussed  by  Dooge  (ibid.,  p.  11). 

Numerous  requests  have  been  received  for  a  description  of  the  detailed 
form  of  two-stage  convolution  as  used  to  produce  nonlinear  response  of  a 
watershed  to  a  time-stream  of  inputs.  Several  requests  were  received 
during  the  Workshop  on  Modeling  Agricultural  Chemicals  in  the  Environment 
held  in  Athens,  Georgia,  in  January  1974.  Since  such  computational  details 
would  not  be  acceptable  for  usual  journal  publication,  this  instructional 
paper  was  prepared  for  distribution. 
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It  is  fortuitous  that  the  ARS-sponsored  Dooge  lectures  are  now  avail¬ 
able  in  published  form  to  provide  the  theoretical  basis  for  convolutional 
models.  Dooge  (ibid.,  p.  23)  gives  the  convolution  of  causal  systems  with 
continuous  data  and  finite  input  as  in  equation  1. 

y(t)  =  /  x(T)h(t-T)dT  (1) 

o 

y(t)  is  the  continuous  time  function  of  output,  h(t-T)  is  the  continuous 
system  response  (impulse)  function,  and  x(T)  is  the  continuous  time  function 
of  input,  t  is  real  time  and  T  is  relative  time  within  the  response  function. 

When  input,  response,  and  output  are  defined  at  discrete  time  points, 
normally  at  uniform  interval,  convolution  changes  from  integration  to 

summation  as  in  equation  2. 

t 

Y(t)  =  I  X(T)  H  (t-T+1)  (2) 

T  -1 

Y(t),  X(T),  and  H(t-  T+l)  are  now  discrete  output,  input  and  response  functions, 
t  and  T  become  discrete  time  variates,  which  are  actually  time  interval  number. 
H(t-T+1)  =  0  for  t-T+1  -  0  for  a  causal  system* 

Equation  2  is  a  matrix-vector  multiplication.  Define  X(T)  as  an  input 


vector  X^,  X2,  X3 ,  X^,  and  X^.  Define  H(t-T+1)  as  the  square  matrix  in 
equation  3. 
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With  these  substitutions,  equation  2  becomes  equation 

4  when  expanded. 
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When  using  two-stage  convolution  to  construct  a  storm  hydrograph 
the  first  stage  is  used  to  construct  watershed  unit  response  functions. 

A  characteristic  function  is  convolved  with  a  state  function.  The  charac¬ 
teristic  function  assumes  the  role  of  input,  the  state  function  is  the 
impulse  function  and  the  watershed  unit  response  function  is  output. 

The  characteristic  function  is  a  time-transformation  of  a  watershed 
potential  runoff  map.  It  is  constructed  as  a  histogram  of  runoff  from 
unit  areas  vs.  time-of-travel  classes.  A  grid  system  can  be  superimposed 
on  the  watershed  to  delineate  unit  areas.  Define  the  characteristic  by  its 

histogram  values  Cl,  C2,  C3 ,  C4,  - ,  CB.  The  time-of-travel  class  width 

is  A  t  .  Assume  that  for  each  interval  At  the  watershed  acts  as  a  linear 
reservoir.  Then  the  state  function  is  equation  5. 


S  ( T) 


=  A  e 


-A(T) 


(5) 


A  is  the  time  constant  of  the  reservoir. 

Snyder  and  Asmussen  (ibid.)  show  that  to  change  this  function  for  discrete 
convolution  the  area  under  the  exponential  function  is  computed  for  each 
At  interval.  That  is,  the  integral  of  equation  5  from  T  =  t  to  T=  t+1 
is  calculated.  The  successive  values  for  S(T)  so  computed  are  given  in 
equations  6. 
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=  S  (2) 

(6-2) 

-2A  _  e-3A 

=  S  (3) 

(6-3) 

-3A  _  e'4A 

=  S  (4) 

(6-4) 

-4A  _  e-5A 

=  S  (5) 

(6-5) 

Generalize  these  values  to  Sl(l),  Sl(2),  Sl(3),  - ,  SI ( 5 )  to  indicate 

these  are  the  values  of  the  state  function  for  the  first  interval  of  time. 

Now  route  each  bar  or  step  of  the  characteristic  function,  as  in  equa¬ 
tions  7,  using  the  form  indicated  in  equation  4.  Only  5  steps  are  given  for 
illustration,  however,  the  form  of  calculation  is  simple  and  easily  extended. 


S1(1)C1 

=  qi (1) 

(7-1) 

S1(2)C1  + 

S1(1)C2 

=  ql (2) 

(7-2) 

S1(3)C1  + 

S1(2)C2  +  S1(1)C3 

=  ql (3) 

(7-3) 

S1(4)C1  + 

S1(3)C2  +  S1(2)C3  + 

S1(1)C4 

=  ql (4) 

(7-4) 

SI (5)C1  + 

SI (4)C2  +  S1(3)C3  + 

S1(2)C4  +  SI (1)C5 

=  qK5) 

(7-5) 

In  equations  7,  we  have  formed  ql(t),  the  watershed  unit  response  function 
for  the  first  interval.  Any  rainfall  excess  generated  in  this  interval  is 
distributed  by  ql(t)  to  form  discharge. 

The  distribution  of  rainfall  excess  is  second  stage  convolution.  For 
this  second  stage, rainfall  excess  is  input,  the  watershed  unit  response 
function  is  impulse,  and  stream  discharge  is  output.  Perform  the  first 
operation  (first  interval  distribution)  of  the  second  stage  convolution 
as  in  equations  8.  Rl  is  rainfall  excess  for  period  1. 


f  - 
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ql(l)Rl  =  Qx 

(8-1) 

ql(2)Rl  -  xQ2 

(8-2) 

ql(3)Rl  =  1Q3 

(8-3) 

ql(4)Rl  =  XQ4 

(8-4) 

ql(5)Rl  =  1Q5 

(8-5) 

In  equations  8,  Q-^  is  calculated  discharge  at  the  end  of  interval 
number  1.  ^Q2  is  calculated  discharge  at  end  of  interval  number  2  from 

rainfall  in  interval  1.  3^3  is  partial  discharge  at  end  of  interval  2, 

and  so  on. 

Consider  the  state  or  routing  function  in  more  detail.  To  get  ql(T) 
in  equations  7,  we  needed  S1(T).  These  values  we  got  from  equations  6, 
by  assuming  we  had  some  value  of  the  reservoir  time  constant,  A.  Suppose 
we  compute  A  as  in  equation  9,  making  it  vary  in  time  but  not  directly  a 
function  of  time. 


A(t)  =  U  +  V  (Qt_±  +  Bt_1)  (9) 

In  this  equation,  Qt_i  is  discharge  at  the  beginning  of  each  discrete 
interval,  t.  Bt_-L  is  antecedent  base  flow,  projected  under  the  storm. 

U  and  V  are  regression-type, or  empirical,  parameters.  At  the  beginning  of 
the  first  interval  storm  discharge  is  zero. 

It  is  well  known  that  watersheds  do  not  act  as  linear  reservoirs. 

When  streams  are  flowing  full,  the  response  to  additional  rainfall  inputs 
is  flashier  than  when  channels  are  nearly  empty.  By  computing  A  as  in 
equation  9,  using  previously  calculated  discharges  as  a  feedback  variable, 
we  can  control  the  flashiness  of  the  watershed,  making  it  act  as  a  nonlinear 
reservoir.  Equation  8-1  gives  the  discharge  at  end  of  interval  1,  which 
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is  obviously 

the  discharge  at 

the  beginning  of  interval  2. 

Specifically 

equations  10 

give  the  values 

of  A  for  the  first  and  second 

intervals . 

A(l) 

=  U  +  V(BC) 

10-1 

A(2) 

=  U  +  VCQj+B-l) 

10-2 

In  retrospect,  A(l)  was  used  in  equations  6  to  get  values  of  S1(T).  This 
reservoir  value  feeds  through  equations  7  and  8  to  give  stream  discharge. 

Now  moving  computationally  forward,  A(2)  will  provide  values  of 
S2(T)  in  equations  6.  The  watershed  unit  response  function  for  interval  2 
is  calculated  in  equations  11. 


S2(1)C1 

= 

q2  (1) 

(11-1) 

S2(2)C1  +  S2(1)C2 

= 

q2(2) 

(11-2) 

S2(3)C1  +  S2 (2) C2  + 

S2(1)C3 

= 

q2(3) 

(11-3) 

S2(4)C1  +  S2(3)C2  + 

S2(2)C3  +  S2(1)C4 

= 

q2(4) 

(11-4) 

S2(5)C1  +  S2(4)C2  + 

S2(3)C3  +  S2(2)C4  + 

S2(1)C5 

q2(5) 

(11-5) 

Proceeding  on, 

the  rainfall  excess 

in  the  second 

interval,  R2 

,  is 

distributed  by  q2(T) 

in  the  second  stage 

convolution. 

Performing 

this , 

and  combining  with  equations  8,  we  get  the  partial  discharges  through 
interval  2,  as  in  equations  12. 


ql(l)Rl 

= 

ii 

(12-1) 

ql(2)Rl  4-  q2(l)R2 

= 

q2 

(12-2) 

ql (3)R1  +  q2(2)R2 

= 

2Q3 

(12-3) 

ql (4)R1  +  q2(3)R2 

= 

2Q4 

(12-4) 

ql (5)R1  +  q2(4)R2 

= 

2Q5 

(12-5) 
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Now  we  can  compute  A(3)  as  in  equation  13. 

A(3)  =  U  +  V(Q2+B2)  (13) 

Equation  6  gives  us  S3(T),  and  we  compute  q3(T)  by  the  same  general 
form  as  equations  7  or  11.  Again,  moving  forward,  we  distribute  the 
rainfall  excess  of  the  third  interval,  R3 ,  and  add  to  what  is  already 
computed,  as  in  equations  14. 


ql(l)Rl 

= 

Q1 

(14-1) 

ql (2)R1  +  q2(l)R2 

= 

Q2 

(14-2) 

ql(3)Rl  +  q2(2)R2  +  q3(l)R3 

= 

q3 

(14-3) 

ql(4)Rl  +  q2 (3) R2  +  q3(2)R3 

= 

3% 

(14-4) 

ql(5)Rl  +  q2(4)R2  +  q3(3)R3 

= 

3^5 

(14-5) 

It  can  be  seen  that  now  A(4) 

can  be 

computed,  response 

functions 

calculated,  and  then  R4  can  be  distributed.  The  process  now  keeps  going 
until  all  rainfall  excess  values  are  distributed. 

The  computations  outlined  above  can  be  performed  by  numerical  methods 
other  than  convolution.  The  reason  for  using  convolution  is  that  it  is 
an  efficient  and  fast  process  in  high-speed  calculation.  A  few  nested 
DO-loops  accomplish  all  the  steps.  Keeping  our  models  short,  fast,  and 
efficient  is  a  fundamental  requirement.  We  must  be  able  to  generate  not 
one  storm  hydrograph,  but  hundreds  or  thousands,  under  many  probabilistic 
sequences  of  future  rainfall  for  many  resource  program  alternatives.  The 
desire  for  extreme  precision  in  calculation  of  one  storm  must  always  be 
tempered  by  the  knowledge  that  there  will  be  probabilistic  variabilities 
among  the  many  storms  of  the  future. 
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A  model  structured  on  convolution,  causes  the  old  but  still  live 
arguments  of  plots  vs.  watersheds  to  assume  proper  perspective.  Consider 
a  plot  within  a  watershed.  This  plot  has  runoff  potential  which  affects 
the  characteristic  function  in  a  watershed.  The  plot  simply  "maps  into" 
the  characteristic  function  histogram,  along  with  all  other  subdrainages 
of  the  watershed.  The  plot  will,  specifically,  affect  some  ordinates  of 
the  characteristic  function.  First  stage  convolution,  equations  7  and  11 
for  example,  shows  how  these  affected  ordinates  would  affect  the  response 
functions.  In  second  stage  convolution,  the  plot  effect  simply  follows 
through  in  the  formation  of  the  discharge  hydrograph. 

If  a  plot  were  changed  so  that  its  runoff  potential  would  change, 
then  a  new  characteristic  function  is  developed.  Beyond  this,  the  plot- 
change  simply  follows  through  the  computations  of  two-stage  convolution 
to  give  the  changed  storm  hydrograph.  While  the  details  of  functional 
formulations  are  open  to  question,  the  basic  computational  perspective 
of  plots  within  watersheds  is  not. 

EPA  is  interested  in  projecting  from  watersheds  to  basins.  The  same 
computational  system  presented  above  applies  just  as  well  in  watersheds 
vs.  basins  as  in  plots  vs.  watersheds.  Again  details  of  formulation  must 
be  worked  out,  but  the  greater  computational  scheme  must  hold.  For 
example,  one  detail  would  be  to  incorporate  the  great  amount  already  known 
about  hydrodynamic  channel  routing. 


APPENDIX 


Given  in  storage: 

Characteristic  function  C(I) ,1=1, NORD 
Rainfall  excess  vector  R(I) ,I=1,IDUR 
Base  flow  at  beginning  of  storm  BO 
Base  flow  under  storm  B(I) ,1=1, NORD 
Number  of  ordinates  NORD 
Number  of  rainfall  intervals  IDUR 
Length  of  interval  DELT  hours 

Compute : 

State  functions  SF(I,J) ,1=1, NORD, J=l, IDUR 
Unit  response  functions  UQ(I,J) ,1=1, NORD, J=l, IDUR 
Storm  discharge  Q(I)J=1,N0RD 
Total  discharge  TQ(I)jI=l  ,N0RD+1 

DO  1000  JB=1 ,N0RD 
DO  1000  JA=1 , IDUR 
1000  UQ (JB , JA) =0 . 0 
A=U+V*B0 
AS  TART =1 . 0 
DO  1001  IT=1 ,N0RD 
AEX=-A*IT*DELT 
AJB=EXP (AEX) 


SF (IT , 1) =ASTART-AJB 


1001 


ASTART=AJB 


DO  1002  KD=l,NORD 
DO  1002  JB=KD ,NORD 
JBI=JB-KD+1 

1002  UQ (JB , 1)=UQ (JB , 1)+SF (JBI , 1) *C (KD) 

DO  1003  JB=1 ,NORD 

1003  Q(JB)=UQ(JB,1)*R(1) 

DO  1004  ID=2 , IDUR 
A=U+V* (Q(ID-1)+B(ID-1) ) 

ASTART=1 . 0 
JBEND=NORD-ID+l 

DO  1005  IT=1 , JBEND 
AEX=-A*IT*DELT 
AJB=EXP (AEX) 

SF (IT , ID) =ASTART-AJB 

1005  ASTART=AJB 

DO  1006  KD=1, JBEND 
DO  1006  JB=KD, JBEND 
JBI=JB-KD+1 

1006  UQ ( JB , ID) =UQ ( JB , ID)+SF (JBI , ID) *C (KD) 
DO  1007  JB=ID ,NORD 

JBI=JB-ID+1 

1007  Q (JB) =Q (JB)+UQ(JBI , ID) *R(ID) 


.  pu  (IeSL)pU 

eoojt  oci 

•ia*TI*A~=X3A 

Jt-KDt-aL=I3L 


1004  CONTINUE 
TQ(1)=B0 
TN=N0RD+1 
DO  1008  JB=2,TN 
1008  TQ (JB  )=Q (JB-1)+B (JB-1) 
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